Višja dinamika

Sistemi z več prostostnimi stopnjami

Dr. Janko Slavič, April 2014

Uporaba fizikalnih koordinat - primer

Remove["Global`*"]

Lagrageva energijska funkcija

L[t_]:=Ek[t]-Ep[t]
Ek[t_]:=(1/2) 6J ((α'[t])^2)+(1/2)J ((β'[t])^2)+(1/2) J ((γ'[t])^2)
Ep[t]:=(1/2)3k ((φ1[t])^2)+(1/2)k ((φ2[t])^2)

Kjer velja:

φ1[t_]:=α[t]-β[t];
φ2[t_]:=β[t]-γ[t]

Uporaba fizikalnih koordinat - gibalne enačbe

Gibalne enačbe

ge1=Simplify[(∂_t) ((∂_(α'[t]))L[t])-(∂_(α[t]))L[t]==0]
ge2=Simplify[(∂_t) ((∂_(β'[t]))L[t])-(∂_(β[t]))L[t]==0]
ge3=Simplify[(∂_t) ((∂_(γ'[t]))L[t])-(∂_(γ[t]))L[t]==0]

Nastavek

α[t_]:=A0 Sin[ω t];β[t_]:=B0 Sin[ω t];γ[t_]:=C0 Sin[ω t]

ge11=ReplaceAll[ge1, Sin[t ω]->1]
ge22=ReplaceAll[ge2, Sin[t ω]->1]
ge33=ReplaceAll[ge3, Sin[t ω]->1]

Uporaba fizikalnih koordinat - harmonski odziv

Zapišemo v matrični obliki

b={A0,B0,C0};
{v,A}=Normal[CoefficientArrays[{ge11,ge22,ge33},b]];

MatrixForm[A]

Poiščimo lastne frekvence

solω0=Solve[Det[A]==0/.ω->ω0,{ω0}]

Pripravimo pravila za zamenjavo lastnih frekvenc

zamenjajω0={ω01->ω0/.solω0[[2,1]],ω02->ω0/.solω0[[4,1]],ω03->ω0/.solω0[[6,1]]}

Uporaba fizikalnih koordinat - lastni vektorji

solAmp1=Solve[(A/.ω->ω01/.zamenjajω0).b==v,{A0,C0}]/.B0->1

solAmp2=Solve[(A/.ω->ω02/.zamenjajω0).b==v,{A0,C0}]/.B0->1

solAmp3=Solve[(A/.ω->ω03/.zamenjajω0).b==v,{A0,C0}]/.B0->1

Uporaba fiziklanih koordinat - lastni vektorji

Določimo lastne vektorje

v1=b/.Flatten[{solAmp1,B0->1}];v2=b/.Flatten[{solAmp2,B0->1}];
v3=b/.Flatten[{solAmp3,B0->1}];

v3

Show[Block[{$DisplayFunction=Identity},
{ListPlot[v1,Joined->True,PlotStyle->Red],
ListPlot[v2,Joined->True,PlotStyle->Green],ListPlot[v3,Joined->True,PlotStyle->Blue]}],Frame->True,PlotRange->All,Axes->False]

Uporaba fizikalnih koordinat - odziv posamezne lastne oblike

odziv1[t_]=v1(C11 Cos[ω01 t]+C12 Sin[ω01 t])
odziv2[t_]=v2(C21 Cos[ω02 t]+C22 Sin[ω02 t])
odziv3[t_]=v3(C31 Cos[ω03 t]+C32 Sin[ω03 t])

Uporaba fizikalnih koordinat - superpozicija odzivov

Sestavljen odziv (vektor)

odziv[t_]:=odziv1[t]+odziv2[t]+odziv3[t]

Dodajanje začetnih pogojev za zasuk (poskusite različne vrednosti za začetne zasuke)

zasuk0={1,1,-2};
solC1=Solve[zasuk0==Chop[odziv[0.]/.zamenjajω0],{C11,C21,C31}]

Dodajanje začetnih pogojev za kotno hitrost (začetne hitrosti so nič)

solC2=Solve[{0,0,0}==Chop[odziv'[0.]/.zamenjajω0],{C12,C22,C32}]

sol=odziv[t]/.Flatten[{solC1,solC2,C12->0,zamenjajω0}/.{k->1,J->1}]

Plot[
{sol[[1]],sol[[2]],sol[[3]]},{t,0,10},PlotStyle->{Red,Green,Blue}]

Glavne ali modalne koordinate - ortogonalnost

MatrixForm[A]

Masna in togostna matrika

M={{6 J,0,0},{0,J,0},{0,0,J}};
K={{3,-3,0},{-3,4,-1},{0,-1,1}}k;

Pokažemo ortogonalnost

Simplify[{v1}.M.Transpose[{v2}]]
Simplify[{v1}.M.Transpose[{v3}]]
Simplify[{v2}.M.Transpose[{v3}]]

Glavne ali modalne koordinate - modalna masa, modalna togost

{{m11}}=Simplify[{v1}.M.Transpose[{v1}]];
{{m22}}=Simplify[{v2}.M.Transpose[{v2}]];
{{m33}}=Simplify[{v3}.M.Transpose[{v3}]];
N[{m11,m22,m33}/.J->1]

Modalne togosti

{{k11}}=Simplify[{v1}.K.Transpose[{v1}]];
{{k22}}=Simplify[{v2}.K.Transpose[{v2}]];
{{k33}}=Simplify[{v3}.K.Transpose[{v3}]];
N[{k11,k22,k33}/.k->1]

Skakalec  - determinanta

Remove["Global`*"]
M={{m,0},{0,m}};
K={{k,-k},{-k,2k}};

Z uporabo determinante

res=Solve[Det[K-(ω^2)M]==0,ω]

rr=N[Chop[ω/.{res/.{k->100 (10^3),m->10}}]]

Skakalec  - lastne vrednosti

Reševanje lastnih vrednosti

M=M/.{k->100 (10^3),m->10};
K=K/.{k->100 (10^3),m->10};

K

B=N[Inverse[M].K]

{values,vectors}=Eigensystem[B];

Take[(Sqrt[Chop[values]])]

Ortogonalnost

ortoM[i_,j_]:={vectors[[i]]}.M.Transpose[{vectors[[j]]}]
ortoK[i_,j_]:={vectors[[i]]}.K.Transpose[{vectors[[j]]}]

ortoM[1,2]

Skakalec  - generalizirana masa in togost

Zaradi vrstičnega zapisa vektorjev, je potrebno vectors transponirati

gM=Chop[vectors.M.Transpose[vectors]]

gK=Chop[vectors.K.Transpose[vectors]]